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Abstract 

Large  temperature  differences  applied  to  thermoelectric  generators  require  that  the  variations  of  all  material  properties  with 
temperature  are  included  in  a  numerical  description  of  their  performance.  A  finite  element  algorithm  is  developed  to  calculate  the 
temperature  field  in  a  thermoelectric  device  and  concomitantly  its  thermoelectric  performance  under  operation  conditions. 
Spatially  varying  the  composition  of  material  or  the  doping  concentration  allows  to  enhance  the  power  output  or  efficiency.  The 
choice  of  the  optimum  concentration  parameter  profiles  is  shown  not  to  be  a  function  of  the  local  temperature  only,  but  to  be 
dependent  on  a  local  criterion  related  to  the  entire  temperature  field.  This  criterion  is  included  in  an  iterative  calculus  to  find 
optimised  concentration  profiles.  It  is  shown  that  the  code  developed  has  advantages  over  previously  published  solutions,  since  it 
can  be  applied  to  continuous  and  discontinuous  material  changes  without  any  assumptions  on  the  mutual  dependence  of  the 
governing  transport  parameters  or  a  need  to  cast  their  temperature  dependence  into  analytical  form.  The  performance  is  shown 
for  some  test  situations  and  compared  to  literature.  ©  1998  Published  by  Elsevier  Science  S.A.  All  rights  reserved. 
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1.  Introduction 

Thermoelectric  generators  are  driven  with  large  tem¬ 
perature  differences  to  convert  the  heat  flow  through 
them  into  electric  energy  with  high  efficiency.  In  typical 
applications  the  difference  across  the  thermoelectric 
material  amounts  to  some  100  K.  A  temperature  inter¬ 
val  as  large  as  this  does  not  usually  allow  consideration 
of  the  material  properties  determining  performance  (i.e. 
Seebeck  coefficient,  thermal  and  electrical  conductivity) 
as  independent  of  temperature.  Because  of  this  fact,  the 
constant  property  approach  commonly  used  in  the  de¬ 
scription  of  thermoelectric  devices,  [1],  is  always  an 
approximation. 

Several  calculation  approaches,  that  include  the  tem¬ 
perature  dependence  of  the  material  properties,  have 
been  published  in  the  literature,  but  there  is  always  a 
need  to  express  the  temperature  dependence  analytically 
or  even  the  mutual  relationships  between  the  material 
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properties  in  order  to  reduce  the  number  of  parameters, 
[3,4,11],  In  particular,  none  of  the  approaches  are  able 
to  deal  with  arbitrary  temperature  dependencies  of  the 
material  properties  which  include  abrupt  property 
changes  as  occurring  in  a  stacked  thermoelectric  ele¬ 
ment  (Fig.  1). 

Here  we  present  an  exact  and  flexible  procedure, 
based  on  a  finite  element  approach.  It  is  easy  to  use  and 
needs  only  tabulated  values  of  material  parameters  as 
input.  There  is  no  physical  restriction  on  their  tempera¬ 
ture  dependence.  Both,  continuous  and  discrete  mate¬ 
rial  transitions  are  correctly  dealt  with.  Even  artificial 
materials  can  be  modelled,  which  offers  the  possibility 
to  deal  with  the  thermoelectric  performance  of  quan¬ 
tum  wells  and  other  new  materials. 

Since  the  material  parameters  usually  depend  on  the 
concentration  of  the  dopants  or  the  composition  of  the 
semiconductor,  it  has  been  suggested  to  exploit  the 
concentrational  dependence  of  one  or  more  of  the 
mentioned,  temperature-dependent,  transport  proper¬ 
ties  to  optimise  the  performance  of  the  whole  device 
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Fig.  1.  A  thermoelectric  generator  is  typically  operated  in  large  temperature  differences.  Graded  or  stacked  devices  improve  the  performance  by 
selecting  the  material  in  accordance  with  the  temperature  field. 


[2,3],  The  basic  idea  is  indeed  simple:  If,  in  a  one-di¬ 
mensional  description,  the  temperature  field  in  a  ther¬ 
moelectric  leg  is  known,  one  can  in  principle  adjust  the 
local  composition  in  order  to  maximise  the  overall 
efficiency  or  output  power  of  the  device  Fig.  1.  How¬ 
ever,  any  alteration  of  the  compositional  profile  will 
change  the  temperature  field  inside  the  leg  in  response, 
which  has  to  be  taken  into  account  in  numerical 
considerations. 

The  maximisation  itself,  i.e.  the  calculation  of  the 
spatial  material  or  compositional  profile,  has  previously 
been  performed  by  either  testing  discrete  materials  at 
spatial  intervals  [9],  or  selecting  compositional  profiles 
[3],  out  of  a  heavily  restricted  variety  with  only  few 
(two  or  three)  free  coefficients,  or  by  applying  standard 
maximisation  procedures  to  the  assumed  functional 
temperature  dependence  of  the  transport  coefficients 
[5], 

Our  practical  solution  deals  with  this  problem  by  an 
iteration  of  the  sequence  ‘calculation  of  the  tempera- 
ture-field/selection  of  material’;  where  a  new  local  crite¬ 
rion  to  select  the  material  composition  is  used. 

It  will  be  shown,  that  the  local  choice  of  material 
always  has  to  consider  the  whole  thermal  environment 
of  each  volume  increment.  In  other  words,  the  local 
contribution  to  the  output  power  is  not  only  a  function 
of  local  temperature  and  composition,  but  depends  on 
the  temperature  field  throughout  the  device.  This  is  in 
contrast  to  known  concepts  that  focus  on  a  ‘local  figure 
of  merit’,  which  may  be  defined  by  Z  =  S2-  a  Ik,  (S, 
Seebeck  coefficient;  a,  electrical  conductivity;  k,  ther¬ 
mal  conductivity).  Here  Z  is  a  material  property.  Max¬ 
imising  Z  at  any  place  has  no  physical  justification  to 
maximise  the  power  output  or  the  efficiency.  As  a  first 
evidence  to  this  point,  one  can  judge  that  Z  is  com¬ 
monly  introduced  in  the  description  of  a  constant  prop¬ 
erty  thermoelectric  generator.  In  a  generator  with 
spatial  variation  of  material  (graded  or  stacked  genera¬ 
tor  as  sketched  in  Fig.  1)  the  situation  is  different:  A 


thin  slice,  where  the  properties  may  be  regarded  to  be 
constant,  is  not  operated  as  a  stand  alone  generator. 
The  electric  current  passing  through  this  single  slice 
arises  from  thermovoltages  produced  by  the  whole 
stack  of  slices  and  is  additionally  determined  by  a  load 
resistor,  which  is  matched  to  the  whole  generator.  The 
influence  of  the  single  slice  on  the  integral  current 
generation  is  only  of  the  order  1  /N  ( N ,  total  number  of 
slices).  Thus  the  local  contribution  has  to  be  a  function 
of  the  whole  thermovoltage  collected  and  the  outer 
load,  which  both  are  integral  quantities. 

The  paper  is  arranged  as  follows:  In  Section  2  the 
finite  element  procedure  for  the  temperature  field  calcu¬ 
lation  is  introduced,  which  is  the  base  to  calculate  the 
performance  of  an  arbitrary  thermoelectric  generator. 
For  practical  applications  it  is  of  great  interest  to 
provide  the  optimum  design  for  graded  or  stacked 
generators  (i.e.  continuous  or  discontinuous  material 
changes).  We  combined  the  calculus  for  the  tempera¬ 
ture  field  with  the  procedure  to  select  the  concentration 
parameter,  Section  3,  and  obtained  suggestions  for 
optimum  device  designs.  Before  presenting  some  results 
of  device  designs  in  the  last  chapter  6,  the  peculiarities 
of  the  inclusion  of  discontinuous  material  changes  into 
the  calculation  are  discussed,  Section  4.  In  Section  5  the 
use  of  this  approach  as  a  design  tool  is  discussed  and 
compared  with  the  literature. 


2.  Calculating  the  temperature  field 

The  temperature  field  inside  a  thermoelectric  genera¬ 
tor  can  be  derived  from  the  known  set  of  coupled 
integral  differential  equations,  [6]: 

0  =  V  •  (kVT)  -  Ti  ■  VS  +  pi 2  =  D([T(x)\),  (1) 
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Table  1 

Terminology  used  throughout  this  article 


Symbol 

Quantity 

Unit 

Fundamental  relation 

X 

Spatial  coordinate 

m 

T 

Temperature 

K 

s 

Flux  of  entropy 

J  K^1  s'1  m-2 

u 

Flux  of  energy 

J  s”1  m-2 

n 

Density  of  particles 

m~3 

j 

Flux  of  particles 

s-1  m-2 

i 

Electric  flux 

C  s-1  m-2 

<f> 

Electrostatic  potential 

j  c-1 

a 

Electric  conductivity 

Q-1  m-1 

P 

Electric  resistivity 

£2  m 

p  =  a~ 1 

K 

Thermal  conductivity 

W  K1  m-1 

S 

Seebeck  coefficient 

V  K-1 

9 

Fleat  flux 

J  s”1  m^2 

q=Ts 

R 

Load  resistance/generator  cross  sectional  a: 

rea  Q  m~2 

R/A 

P 

Energetic  potential  of  the  charge  carriers 

J 

(For  the  meaning  of  the  symbols  refer  to  Table  1.) 
These  equations  state  the  conservation  of  energy  and 
number  of  charge  carriers  as  well  as  linear  transport 
obeying  the  Onsager  relations. 

To  solve  Eq.  (1),  we  apply  an  iterative  approach  to 
take  the  integral  character  of  Eq.  (2)  into  account.  This 
allows  to  treat  the  coefficients  of  Eq.  (1)  just  as  func¬ 
tions  of  temperature,  but  independent  of  the  whole 
temperature  field.  We  restrict  the  solution  to  one-di¬ 
mensional  fields,  i.e.  T  =  T(x),  etc.  because  we  focus  on 
large  temperature  differences  along  the  axis  of  a  ther¬ 
mogenerator.  If  there  are  transverse  differences,  they 
are  due  to  inhomgeneous  heat  flux  at  the  hot  or  cold 
junction  of  the  material  or  due  to  radiation  losses. 
Typical  devices  are  designed  to  have  isothermal  junc¬ 
tions,  therefore  any  transverse  temperature  variation 
will  be  much  smaller  than  the  axial  one. 


Fig.  2.  Temperature  profile  in  the  n-leg  of  a  thermoelectric  generator. 
The  temperature  dependent  material  properties  were  introduced  by 
Sherman  et  al.  [4],  Within  the  precision  of  the  results  available  from 
Ref.  [4]  the  results  match  perfectly.  The  inset  shows  the  remaining 
difference  between  both  profiles. 


Thermoelectric  devices  are  generally  assembled  form¬ 
ing  pairs  of  p  and  n-type  material,  the  p-  and  n-leg  with 
positive  and  negative  Seebeck  coefficients,  respectively, 
Fig.  1.  Often  such  a  couple  is  treated  as  a  unit  [4], 
where  both,  the  p  and  the  n-leg,  are  included  in  the 
calculus.  But  it  is  sufficient  for  a  comprehensive  charac¬ 
terisation  and  optimisation,  to  treat  the  legs  individu¬ 
ally.  It  can  be  shown,  that  optimised  conditions  require 
a  matched  electric  current  density.  When  combining  the 
two  legs  to  the  desired  thermocouple,  the  respective 
current  densities  in  the  two  materials  have  to  be  con¬ 
served  by  adjusting  the  area  ratio  of  the  legs  and 
combining  the  two  calculated  load  resistances. 

Under  these  conditions  a  one  dimensional  tempera¬ 
ture  field  is  calculated  using  a  direct  solution  according 
to  Galerkin’s  method  of  the  finite  element  algorithm 
(e.g.  [7,8]).  We  used  equidistant,  linear  element  func¬ 
tions  f  (x)  along  the  material.  The  temperature  field 
thus  has  the  following  representation: 

7'M  =  X  Ttf(x).  (3) 

;,xe[0,i] 


Table  2 

Design  of  a  stacked  generator.  3  (2)  hypothetical,  constant  property 
materials  are  combined  to  form  the  p  (n-)leg,  resp.,  of  a  generator 


Material 

designation 

Temperature 
range  [K] 

Relative  length 
of  the  segment 

Reference  [11] 

This  work 

PI 

880... 720 

0.551 

0.575 

P2 

720... 420 

0.291 

0.280 

P3 

420...  300 

0.158 

0.145 

N1 

880... 420 

0.772 

0.795 

N2 

420...  300 

0.228 

0.205 

The  results  of  Ref.  [1 1],  who  have  proposed  the  material  properties, 
are  close  to  this  approach. 
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Fig.  3.  Calculating  the  optimum  concentration  profile  in  a  twofold  stacked  generator.  Provided  all  other  properties  are  identical,  the  substance 
with  maximum  S\  .is.  selected.  The  two  available  S(T)  are  shown  in  the  left  part.  Temperature,  T(x)\  temperature  gradient,  VT(x)  and  material 
parameter  c(x)s{0,  1}  are  collected  at  the  right. 


The  profiles  of  the  material  parameters  are  repre¬ 
sented  in  the  same  way.  E.g. 

K(x)  =  X  Ki  ~  k(T  =  T(x  =  iL/N)). 

i,xe[0,L] 

According  to  the  Galerkin  method  [7,8],  a  set  of 
equations,  that  allows  the  determination  of  the  {7V},  is 
obtained  from 

j7  Ux)D({T(x)])dx  =  0,  i=l...N  (4) 

With  the  known  functions  f  (x)  the  integration  of 
Eq.  (1)  is  straightforward. 

The  transport  coefficients  S,  a  and  k  are  evaluated, 
considering  the  temperature  T  and  the  concentration  c 
as  independent  variables  and  using  the  values  of  the 
previous  step  of  the  iterative  calculation  of  T(x).  All 
derivatives  are  expressed  via  the  T  and  c  fields: 

e.g,V«.^V7-  +  |v,  (5) 


Thus  the  material  parameters  S,  a  and  k  and  their 
derivatives  can  easily  be  interpolated  from  tables 
S(T,  c ),  a(T,  c ),  and  k(T,  c),  that  are  given  to  this 
routine  as  input  files.  We  consider  this  a  simpler  way 
than  eliminating  the  temperature  or  concentration 
parameter  using  analytical  relations  [3,5], 

Since  Eq.  (1)  is  of  second  order,  a  quadratic  set  of 
relations  between  the  coefficients  {7}}  is  obtained  from 
Eq.  (4).  For  one-dimensional  T{x),  as  considered  here, 
one  single  equation  for  Ti+1  needs  to  be  solved,  once 
Tt_x  and  Tt  are  known,  or  equivalently  T  and  VT  at 
one  boundary.  The  corresponding  quadratic  equation 
and  the  appearing  coefficients  are  given  below  for  the 
general  case  of  temperature  and  concentration  depen¬ 
dent  S(T,  c ),  <r(7;  c)  and  k(T,  c). 

0=Tj+l-hl  +  Ti+l(h2  +  h3Ti) 

+  (h4  r?  +  h5  Tt  +  h6  Tt  Tf_  J  +  h7  Tj_ ,  +  hs  Tf_ ! 

+  h0),  (6) 


1  (dK\  _  <5/c I  \  i  f  <3S|  051  \ 

1  65x2^1,.+!  +  2  dT\ J  +  iMxVarl,. +  02V  J 

,  Ki  i  f8S\  0S|  05 1  0S| 

2  dx2+ 12dxVdc|,C!+1  0c|,-C,+ 0c|,-+iC,  +  1  dc|i+1C' 

1  /  „  Sk\  _  0jc|  Sk\  0k I 

+  — jl  ~2^\  ci  +  2^\  ci+ c<  +  —  ct 

OQX  \  uC  OC  |z-  OC  1 

i  05 1  1  /  0k|  Sk\  \ 

3  =  6d.v  8T\t  ~  3d.v7V  r?7  , +  df\ i+1f 


where 
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jvfkl  5k  I  5k  I  \  _i_/5S|  _dS_\  Y 

Mx\8T\i+ ,  +  4  5T|f  +  dl\,_  J  +  12d.v\f;7'|t  ,  dT\i+  J’ 

-  2  K{  i(  -  5S/5c|,_  i  c,_ !  4  ®Sy5c|,_*  %  -  3  8S/dc\,  ct_  j  +  3  dS/dc\,  ci+ 1  -  dS/dc\i+ ,  %  +  5S/5c|/+ *  c,  +  0 

dx2  +  1 2d  v 

(  —  dic/dc c,_!  +  5/c/5c|f_  x  c,-  —  25k /5c |,-  ct_  x  +  48k  /dc^Cf  —  28k  /8c\icl+x  +  c?K/ocj;  +  1  cf—  5k/5c|,-+1c,-+1) 
+  6d?  ’ 

—  1  /5k|  5k|\  /  5S| 

6d*  Vr|,_ ,  +  2dT\  J6dx  571/ 

1  /5k  I  5k  I  \  -jj  /5S|  55|  \ 

6d.vV7’|,  ,  +  2df\ J  ~  +  cv{  J* 

1  i  (8S\  5S|  55|  8S I  \ 

dx^^ Ki +  12dx(¥|,._  1  +  Jc[_  ,Ci ~  to[Ct- 1  +  Jc[Ci) 

1  /5k|  5k I  _5k|  _5k|  \ 

+ |,_ . c- '  ~  4- . c,+ ^ *4  ■ c- : 1  _  r 


( c  concentration  parameter,  d x  =  L/N) 

With  the  aid  of  this  algorithm  it  is  possible  to 
determine  the  temperature  field  inside  any  given  genera¬ 
tor  leg  under  operation,  if  S,  k,  p  are  known  (i.e. 
tabulated)  functions  of  temperature,  and  concentration 
if  desired. 

Once  having  obtained  the  temperature  field  for  a 
given  condition,  i.e.  material  distribution  and  electric 
load,  it  is  straightforward  to  determine  all  generator 
characteristics.  For  example,  the  internal  resistance  is 
calculated  by  integrating  the  temperature  dependent 
resistivity  along  the  known  temperature  profile;  the 
specific  output  power  is  given  by  the  electrical  current 
and  the  load  resistance;  and  to  receive  the  efficiency  one 
additionally  calculates  the  heat  flow  input  by  evaluating 
the  temperature  gradient  at  the  hot  side.  For  the  sake 
of  simplicity  an  average  value  for  the  gradient  and  the 
temperature  are  taken  as  boundary  conditions  in  the 
first  step.  The  gradient  is  varied  until  the  cold  junction 
temperature  is  matched.  Since  all  thermoelectric  effects 
are  in  the  order  of  10%  (efficiency)  or  less,  the  nor¬ 
malised  difference  between  the  temperature  field  with 
and  without  thermoelectric  effects  (i  =  0)  is  expected  to 
be  of  the  same  order  of  magnitude.  Therefore  a  starting 
value  of  i  =  0  is  chosen,  to  calculate  a  first  approxima¬ 
tion  for  T(x).  This  temperature  field  is  used  to  calculate 
the  next  order  value  for  i.  The  procedure  is  iterated 
until  the  electrical  current  as  well  as  the  temperature 
field  have  stabilised.  We  have  not  experienced  any 
trouble  concerning  convergence. 

We  implemented  the  calculus  on  a  standard  PC  using 
C++  programming.  The  demands  on  computational 
resources  for  ~  100  nodes  are  still  low.  An  example 
calculation  will  be  given  in  Section  6. 

The  additional  concentration  parameter  so  far  only 
serves  to  introduce  different  temperature  dependencies 
at  different  places  resembling  a  certain,  fixed  material. 


In  the  next  section  we  will  show  how  to  select  that 
particular  material  that  maximises  the  electrical  power 
output. 


3.  Spatial  material  selection 


To  maximise  the  power  output  or  efficiency  of  a 
thermoelectric  generator,  one  has  to  make  a  choice 
between  different  materials  at  each  place.  The  efficiency 
of  a  device  is  defined  as 


Pe  i 


(7) 


(Pcb  electric  power  output;  Qin,  total  amount  of  energy/ 
heat  flow  supplied  to  the  generator  at  the  hot  junction.) 
Consider  a  stacked  generator  of  M  different  materials. 
If  there  are  N  slices,  the  optimum  design  has  to  be 
selected  among  MN  different  candidates.  Such  proce¬ 
dures  have  been  employed  by  Anatychuk  et  al.  [9],  but 
for  large  M  and  N,  approaching  continuous  grading, 
this  becomes  intractable  and  a  well  founded  criterion 
for  material  choice  is  needed. 

We  only  consider  the  power  output  here,  because  this 
quantity  is  to  be  maximised  in  most  applications  of 
generators. 

For  an  arbitrary  generator,  the  power  is  given  by  the 
difference  between  outgoing  and  incoming  electric  en¬ 
ergy  flux,  pj.  (The  complete  flux  of  internal  energy  is: 
u=Ts  +  pj) 


Pe  j 


L 


pj-dAi 


pj-dA 

div(pj)dV. 


pj-dA0 


(8) 
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To  evaluate  Eq.  (8),  the  conservation  of  charge  carri¬ 
ers,  n  =  0<s>div/'  =  0,  and  the  representation  of  the  parti¬ 
cle  flu x,  j  =  a/e2(SVT  —  eVp)  [6],  are  used.  As  in  the 
previous  Section,  the  generator  is  considered  to  be  one 
dimensional  with  a  cross  sectional  area  A. 


Pe  i  =  |  +  gSVr)d  v  (9) 

=  A$  dxjejSVT- e2j2p) 

=Aej  p  dT  S  —  Ae2j2  J  d x  p.  (10) 

This  shows  that  the  total  electric  power  output  Pe, 
can  be  written  as  a  volume  integral  over  a  local  contri¬ 
bution  7reI.  Therefore,  the  maximum  power  condition  is 
satisfied,  if  the  integrand  is  maximised  in  every  place. 
This  is  possible  only,  if  the  particle  flux  j  is  known.  It 
can  be  obtained  from  basic  thermodynamics:  The  po¬ 
tential  ji  is  related  to  the  electrostatic  potential  </>  by 
A ^  —  eAcj).  In  an  external  load  circuit  (with  resis¬ 
tance  R)  the  current  is  given  by  I=A</>/R=  —  ejA . 
Thus  An  =  e2jAR.  On  the  other  hand 


/*  =  J  dxV/u=  —  e2j^  d xp  —  ej  d x  SVT. 
j1  dx  SVT 


(ii) 


To  achieve  maximum  power  output  of  a  power  sup¬ 
ply,  which  exhibits  constant  EMF,  independently  of  the 
current  flow,  the  external  load  resistance  R  must  match 
the  requirement  R  =  jo  dx  p.  This  holds  true  for  fixed 
hot  and  cold  side  temperature.  The  right  hand  side  of 
Eq.  (11)  is  thus  determined.  With  the  j  value  obtained, 
the  choice  of  maximum  nA  at  every  position  maximises 
the  output  power. 

7rel  contains  the  integral  quantity  j,  that  characterises 
the  environment,  not  only  the  material.  This  is  indeed  a 
new  fact  in  comparison  with  selection  rules  that  are 
based  on  a  ‘local  figure  of  merit’,  Z  =  ( S2<t)/k ,  employ¬ 
ing  (temperature  dependent)  local  properties. 

Practically,  it  can  be  checked  in  simple  cases  that  a 
material  selection  based  on  Eq.  (10)  does  reproduce  the 
expected  results:  Consider  two  materials  with  equal  and 
constant  thermal  and  electric  conductivity.  The  Seebeck 
coefficient  of  one  species  may  rise  with  temperature 
stronger  than  the  other.  They  are  at  the  same  level  at  a 
certain  temperature  Te.  The  first  term  in  Eq.  (10)  is 
maximized,  by  choosing  the  material  with  maximum  |S|. 
The  second  term  doesn’t  change  with  the  material 
selection,  if  p(T)  is  constant.  Maximising  nel  according 
to  Eq.  (10)  is  to  select  the  material  with  maximum  |S|, 


i.e.  an  interface  between  both  substances  is  thus  to  be 
placed  at  temperature  Te.  As  soon  as  p  actually  de¬ 
pends  on  T,  both  terms  in  Eq.  (10)  have  to  be  consid¬ 
ered.  Then,  the  developed  numerical  method  is 
necessary  for  fixing  the  material  interface  temperature. 

We  implemented  a  computer  code  to  select  optimum 
materials  in  a  given  temperature  field.  This  module  was 
combined  with  the  algorithm  to  determine  T(x).  In  this 
way  a  software  tool  has  been  designed,  that  gives  the 
layout  of  one  thermogenerator  leg  to  produce  maxi¬ 
mum  electric  power  at  fixed  thermal  boundary  condi- 
tions.The  algorithm  follows  the  sequence  below: 

1.  Choose  an  initial  profile  of  the  concentration 
parameter  and  select  the  thermal  boundary  condi¬ 
tions,  Th  and  Tc. 

2.  Choose  external  resistance. 

3.  Calculate  the  temperature  field. 

4.  Using  the  temperature  field  of  3.  and  the  overall 
current  related  to  it,  select  the  optimum  concentra¬ 
tion  parameter  at  all  places. 

5.  Check  whether  the  concentration  and/or  tempera¬ 
ture  have  changed  significantly.  If  not:  quit,  else 
continue  at  2. 

We  experienced  only  two  or  three  loops  to  establish 
a  final  concentration  field  of  three  stacked  materials  in 
a  stack  of  some  hundred  slices.  Altogether  this  allows 
to  make  detailed  suggestions  on  the  design  of  graded  or 
stacked  thermoelectric  generators  with  comparably  lit¬ 
tle  computational  effort. 


4.  Continuous  and  discontinuous  material  changes 

So  far  the  concentration  parameter  c  has  been  treated 
as  a  floating  real  number.  This  is  adequate  for  prob¬ 
lems,  where  c  describes  a  dopant  concentration  as  may 
be  in  FeSi2:Co  (see  example  4  in  Section  6),  or  the 
composition  of  a  matrix  like  in  Ge^Si,  _x.  Here  one 
intends  to  create  a  smooth  transition  or  concentration 
profile. 

The  situation  is  different  in  stacked  generators.  Pro¬ 
vided  no  complications  arise  at  the  junction,  it  is  possi¬ 
ble  to  combine  a  high  temperature  material,  e.g. 
Ge03Si07,  with  medium  temperature  material,  e.g. 
PbTe,  to  increase  the  accessible  temperature  interval.  In 
this  case,  one  asks  for  the  best  position  of  the  interface. 
Of  course,  one  imagines  a  stepwise  material  transition 
for  which  the  basic  representation  of  c(x),  as  given  in 
Eq.  (3),  is  not  appropriate.  We  have  chosen  another 
possible  form  of  the  concentration  profile  for  this  case. 
Instead  of  the  linear  element  functions  in  Eq.  (3), 
constant  elements  are  used: 


c(x).=  Y,  cigi(x)- 

i,xe[0,L] 


where  gt(x)  = 


0,  if  xt[i.L/N,(i+  1  )L/N[, 
1,  if  xe[iL/N,(i+\)L/N[. 


(12) 
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( L ,  overall  length  of  the  leg,  A,  total  number  of 
elements). 

This  alters  the  integrals  of  Eq.  (4).  They  have  to  be 
re-evaluated  with  g,(x)  instead  of  fix),  and  the  corre¬ 
sponding  Eq.  (6)  changes.  But  again  a  quadratic  equa¬ 
tion  for  the  determination  of  7) ,  is  received,  that  has 
the  following  form: 

a  ■  Tj+ 1  +  b-  Ti+ ,  +  c  —  0, 

where  a,  b,  c  are  still  products  of  material  properties, 
evaluated  employing  the  T  field  of  the  previous  step, 
and  powers  of  the  temperature  coefficients  T,  and  TiA. 
Thus  the  subsequent  determination  of  the  temperature 
field  is  not  changed. 

The  selection  procedure  to  determine  the  ‘best’  mate¬ 
rial  in  each  place  is  not  affected  by  the  choice  of 
continuous  or  discontinuous  material  changes.  Using 
linear  element  functions  in  the  representation  of  c(x ), 
the  material  properties  will  vary  linearly.  So  no  extreme 
can  occur  between  two  successive  tabulated  values.  The 
maximum  value  of  nel  in  Eq.  (10)  will  always  be  reached 
at  one  of  the  filed  values  of  the  concentration  parame¬ 
ter  c.  Obviously  the  same  holds  true  in  case  of  a 
stepwise  constant  c(x),  using  Eq.  (12). 

Even  though  we  allow  for  stepwise  material  changes, 
the  optimum  position  of  the  interface  is  only  deter¬ 
mined  to  a  resolution  of  2 L/N,  not  L/N  as  one  might 
have  expected.  The  reason  for  this  is  to  be  seen  in  the 
spatial  dependency  of  the  variations  in  S  and  k:  Con¬ 
sider  a  stepwise  change  in  k  at  a  given  temperature 
Tstep.  Unless  a  node  accidentally  happens  to  fall  to  a 
position  where  T  =  Tstep,  dx/dT  =  0  for  all  nodes. 
Hence  there  is  no  contribution  from  the  terms  regard¬ 
ing  derivatives  of  S  or  k;  in  Eq.  (6)  at  all  nodes.  This 
problem  has  to  be  overcome  by  further  manipulation  to 
the  temperature  Ti+1  calculated  from  Eq.  (6).  The  local 
conservation  of  heat  flux  has  to  be  checked  explicitly. 
In  case  the  incoming  heat  flux,  —  k  (Tt  —  7)_  ,)/dx, 
does  not  match  the  outgoing  flux  —  k  (7)  —  Ti_l)/dx, 
one  can  conclude,  that  a  stepwise  change  of  k  or  S  has 
taken  place.  The  corresponding  value  of  Tt+i  is  ad¬ 
justed  appropriately.  Because  two  neighbouring  tem¬ 
perature  differences  are  considered,  the  spatial 
resolution  for  a  stepwise  change  of  conductivity  is 
decreased  to  2dx. 


5.  Assessment  as  design  tool 

5.1.  Flexibility  of  approach 

The  combination  of  the  finite  element  calculation  of 
the  temperature  field  under  operating  conditions  and 
the  compatible  choice  of  the  optimum  material,  max¬ 
imising  the  local  power  contribution  nA,  gives  a  flexible 


and  easy  to  use  numerical  tool.  It  allows  to  assess  the 
performance  of  existing  devices  and  may  be  used  as  a 
design  tool  when  looking  for  optimised  dopant  profiles 
in  continuously  graded  devices  or  the  best  interface 
position  of  stacked  generator  legs,  incorporating  any 
material  properties  desired. 

This  tool  may  be  compared  to  a  different  one  sug¬ 
gested  by  Mahan  [3]:  Instead  of  the  representation  of  S, 
a,  and  k  as  functions  of  T  and  c,  he  eliminates  the 
concentration  parameter  and  chooses  a  as  independent, 
S  =  S(T,  <t),  k  =  k(T,  a),  by  using  the  Price  relations, 
[10].  After  applying  his  procedure  the  back  transforma¬ 
tion  a(x)  -» c(x)  has  to  be  done  to  receive  a  concentra¬ 
tion  profile,  if  desired  as  design  parameter.  The  calculus 
proceeds  by  expanding  a(x )  as  series  of  Legendre  poly¬ 
nomials.  Similar  to  the  approach  presented  here,  first  a 
temperature  field  for  fixed  a(x )  is  derived  by  variation 
of  the  temperature  gradient  at  one  boundary  to  match 
the  other  boundary  temperature. 

As  in  this  work,  the  determination  of  T(x)  is  part  of 
an  outer  loop  to  find  the  optimum  a  profile.  Since  in 
Ref.  [3]  efficiency,  not  power,  is  maximised,  the  integral 
value  of  the  efficiency  is  directly  used  to  decide  whether 
the  a  profile  has  converged  or  not.  No  local  physical 
criterion  is  used,  see  Section  3. 

5.2.  Numerical  effort 

Employing  our  approach,  the  determination  of  the 
interface  position  in  a  stacked  generator  is  possible  with 
a  relative  precision  of  2/A,  A  number  of  nodes,  see 
Section  4.  Decreasing  the  uncertainty  thus  results  in  an 
increased  numerical  effort  growing  linearly  with  A.  This 
holds  true  for  both,  continuous  and  stepwise  material 
changes,  since  only  the  representation  of  c(x)  is 
changed,  not  the  structure  or  number  of  equations  to 
be  solved. 

In  Ref.  [3]  it  is  claimed  that  only  two  Legendre  terms 
are  required  for  a  fair  representation  of  a(x).  In  this 
case  continuous  grading  is  considered  and  the  computa¬ 
tional  demands  should  be  extremely  low:  For  each 
determination  of  a  T  field  only  the  algebraic  relations 
S(T,  a),  k(T,  a)  have  to  be  evaluated  as  well  as  a 
coupled  set  of  two  equations  for  the  two  Legendre 
coefficients. 

The  representation  of  a  discrete  junction  is  more 
complicated:  To  achieve  the  same  resolution  as  with  the 
finite-element  approach  employing  A  nodes,  the  num¬ 
ber  of  Legendre  terms  required  is  of  order  A.  Solving 
for  the  Legendre  coefficients  results  in  N  equations  with 
roughly  A  terms  each.  Thus  the  numerical  effort  in¬ 
creases  with  the  order  of  A2.  Questions  how  to  derive 
the  individual  coefficients  from  these  equations  are  not 
even  considered  in  this  estimate. 
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Fig.  4.  Recalculation  (subscript  calc)  of  the  device  performance  of  a  FeSi2  generator  analysed  in  Refs  [12,13]  (subscript  exp).  The  simulation 
neglects  the  finite  electric  resistance  of  the  hot  junction,  thus  the  calculated  values  of  the  internal  resistance  are  lower  than  the  experimental  ones. 


5.3.  Describing  different  materials 

Another  interesting  difference  between  Ref.  [3]  and 
this  approach  results  from  the  usage  of  the  Price  rela¬ 
tion  [10],  in  Mahan’s  paper.  The  empirical  relations 
S  =  S(T,  a)  and  k  =  k(T,  a),  will  allow  to  cover  a 
variety  of  materials  and  dopants.  They  use  two  free 
coefficients,  the  energy  gap  and  a  scattering  parameter, 
that  may  be  chosen  for  each  material.  But  still  there  is 
no  physical  need  for  these  relations  to  be  generally 
valid.  Forthcoming  materials,  e.g.  quantum  well  materi¬ 
als  for  S  enhancement,  should  require  additional  free 
parameters  to  find  equivalent  relations.  Even  classical 
counter  doped  material  will  not  fit  into  the  Price 
scheme  as  soon  as  the  temperature  allows  both  dopants 
to  be  significantly  activated  (since  there  is  just  one 
parameter  for  the  energy  gap). 

Regarding  this,  the  use  of  Price  relations  or  similar 
ones  for  different  material  classes  will  always  pose 
limitations  on  the  ‘allowed’  a  interval.  As  soon  as  the 
calculus  of  Ref.  [3]  suggests  values  beyond  this  range, 
the  output  is  no  longer  useful  as  information  for  device 
designers.  We  consider  our  approach  superior  in  this 
point:  The  input  data  are  filed  tables  of  S,  a,  and  k 
versus  T,  c.  The  material  choice  is  restricted  to  the  filed 
values  and  can  be  regarded  as  a  choice  between  exist¬ 
ing,  tested  and  desired  materials.  The  treatment  of  new, 
non-classical  materials  is  easily  included,  once  the  de¬ 
pendencies  S(T,  c ),  a{T,  c),  and  k(T,  c )  are  known. 


6.  Results 

It  is  difficult  to  obtain  reliable  data  on  the  perfor¬ 
mance  of  thermoelectric  materials  for  different  ‘concen¬ 


tration  parameters’,  like  dopant  concentrations  or 
matrix  composition.  So  at  this  stage  we  restrict 
ourselves  to  four  cases:  (1)  Calculation  of  the  tempera¬ 
ture  field  for  one  single  material.  (2)  Determination  of 
the  temperature  field  in  a  stacked  generator.  The  differ¬ 
ent  materials  are  characterised  in  adjacent  temperature 
intervals.  The  temperatures  for  the  material  changes  are 
thus  fixed  and  this  allows  to  treat  this  configuration  as 
one  material  with  stepwise  changing  properties,  skip¬ 
ping  the  selection  procedure.  (3)  Selection  of  one  out  of 
two  materials  with  maximum  Seebeck  coefficient.  (4) 
Simulation  of  experimentally  derived  performance. 

1.  First  we  prove  that  this  approach  allows  to 
calculate  reliable  temperature  fields.  Sherman  et  al.  [4], 
have  established  a  numerical  procedure  to  determine 
the  performance  of  thermoelectric  devices.  We  have 
chosen  their  example  II  for  power  generators  for  a 
comparison.  The  temperature  profiles  are  indeed  very 
close  to  each  other.  Subsequently  all  other  performance 
data  coincide.  We  could  obtain  no  original  numerical 
data  sets  from  their  work,  so  we  digitised  their  pub¬ 
lished  diagrams.  Nevertheless  a  deviation  of  <  1%  is 
found  for  the  local  temperature.  Within  the  limits  of  the 
possible  precision  we  consider  this  conformity  perfect. 
Our  results  for  the  n-leg  are  compared  to  the  data  of 
Ref.  [4]  in  Fig.  2. 

2.  In  the  work  of  Swanson  et  al.  [1 1]  the  design  for 
a  stacked  generator  is  discussed.  They  assume  fixed 
temperatures  for  the  material  changes  and  proceed  as 
follows:  The  interface  positions  between  successive  ma¬ 
terials  are  determined  by  considering  the  net  heat  fluxes 
and  the  electric  output  power.  Concerning  conduction 
of  heat,  a  linear  temperature  profile  is  assumed  in  each 
segment,  neglecting  distortions  due  to  Joule  and  Thom¬ 
son  heats  as  well  as  temperature  dependence  of  thermal 
conductivity.  We  found  the  interface  positions  directly 
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from  the  temperature  profiles.  The  stepwise  changes  of 
material  properties  at  the  given  interface  temperatures 
cause  the  temperature  gradient  to  be  discontinuous. 
The  interface  positions  are  located  at  these  discontinu¬ 
ities.  Since  here  the  fully  non-linear  temperature  profile 
is  included,  slight  deviations  between  Ref.  [11]  and  our 
work  are  expected.  A  fair  coincidence  of  both  ap¬ 
proaches  is  summarised  in  Table  2. 

3.  The  selection  of  the  ‘best’  material  is  demon¬ 
strated  in  a  simple  situation.  A  stacked  generator  leg  of 
two  materials  is  considered.  Both  species  have  identical, 
constant  electrical  and  thermal  conductivity.  The  varia¬ 
tion  of  the  Seebeck  coefficient  with  temperature  is 
shown  in  the  left  diagram  of  Fig.  3.  The  concentration 
parameter  changes  abruptly  from  one  material  to  the 
other  at  a  certain  position.  The  corresponding  interface 
temperature  is  depicted  in  Fig.  3.  One  finds  that  the 
procedure  locates  the  interface  exactly  at  the  tempera¬ 
ture  where  both  materials  have  the  same  Seebeck  coeffi¬ 
cient.  The  material  with  the  maximum  magnitude  of  S 
is  selected,  as  one  has  expected. 

4.  The  group  of  Birkholz  has  published  detailed 
experimental  results  on  the  performance  of  their  gener¬ 
ators  [12,13],  In  these  works,  a  cylindrical  FeSi2-genera- 
tor  is  analysed,  p  and  n-leg  are  connected  by  a  bridge  of 
sintered  FeSi2  of  2.5  mm  thickness.  The  individual  legs 
have  a  length  of  10.5  mm  each.  In  our  calculation  we 
assumed  an  ideal  hot  junction  showing  no  thermal  and 
electrical  resistivity.  The  material  properties  were  only 
documented  for  the  n-material,  FeSi2:Co  [14],  FeSi2:Al 
similar  to  the  p-material  used  in  Refs  [12,13]  was 
analysed  by  Hesse  [15],  Since  the  measurement  of  ther¬ 
moelectric  properties  is  difficult  in  any  case  and  sensi¬ 
tive  to  the  very  material  under  consideration,  the 
agreement  between  experimental  and  calculated  perfor¬ 
mance  data  is  very  satisfactory.  Differences  in  the  resis¬ 
tance  can  be  ascribed  to  the  electrical  junction 
resistance,  which  is  neglected  in  our  calculation  Fig.  4. 

These  examples  clearly  show  that  the  calculation  of 
the  temperature  field  performs  well.  The  subsequent 
choice  of  optimum  concentration  profiles  is  harder  to 
judge  [16],  But  the  simple  tests  presented  here,  repro¬ 
duce  classical  results  exactly  and  the  recalculation  of 
available  studies  on  stacked  generators  gives  satisfying 


agreement.  Though  not  using  the  material  selection 
procedure,  the  simulation  of  experimental  data  on 
whole  generators  shows  the  flexibility  and  reliability  of 
this  calculus. 
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